capture log close 
log using ptverification.log, replace

***************************************************
***    Are coups good for democracy 		***
*** 						***
*** 	author: jgw				***
***	origin date: 8.12.15			***
***	last update: 8.12.15			***
***						***
***	using files: 				***
***		thyne_powell_FPA2013_repdata.dta***
***		GWF_AllPoliticalRegimes.dta	***
***						***
***************************************************

*****************************************************************************************
*											*
* \copyright   2015 Joseph Wright \\							*
*											*
* This research is funded by the National Science Foundation BCS-0904463		*
*											*
* This program is free software: you can redistribute it and/or modify			*
*    it under the terms of the GNU General Public License as published by		*
*    the Free Software Foundation, version 3. \\					*
*											*
* This program is distributed in the hope that it will be useful,			*
*   but WITHOUT ANY WARRANTY; without even the implied warranty of			*
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the			*
*   GNU General Public License for more details.\\					*
*											*
* See  http://www.gnu.org/licenses/ for a copy of the GNU General Public License. 	*    
* 											*
*****************************************************************************************

******************************
** Powell-Thyne replication **
******************************

set matsize 800
cd "C:\Users\jwright\Documents\My Dropbox\Research\Coups\DFGW\Data files\annual data analysis"
use thyne_powell_FPA2013_repdata, clear
qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3, nolog
sum year if e(sample)
recode ccode (679 = 678) (364=365) (529=530) (769=770) (818=816) /* Yemen, Russia, Ethiopia, Pakistan, Vietnam */
rename ccode cowcode
sort cow year
merge cow year using GWF_AllPoliticalRegimes
tab _merge 

egen caseid = group(gwf_case)
gen mil = gwf_military==1 if gwf_mil~=. | dem~=.  /* cover both GWF and PT samples */
recode mil (0=1) if cow==950 & ((year>=1998 & year<=1992) | year==2001 | (year>=2007 & year<=2008))   /* Fiji military dictatorships */
tsset cow year
gen lagmil = l.mil
qui:logit dem coup_a prev_dem colbrit ind cold gdppc chgdppc time time2 time3, nolog
gen s = e(sample)  /* estimating sample */
egen count = count(cow) if s==1 & gwf_mil==., by(cow)
tab cow if s==1 & gwf_mil==. & count>1  /* countries in estimating sample not in GWF */

	* verification model *
	qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3, nolog
	sum year mil gwf_mil if e(sample)
	lincom coup_s
	est store rep1
	* add military regime *
	qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 mil, nolog
	sum year mil gwf_mil if e(sample)
	lincom coup_s
	est store rep2
	* add lagged military regime *
	qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 lagmil, nolog
	lincom coup_s
	est store rep2a
		tab coup_s mil if s==1, col
		tab dem mil if s==1, col
		tab dem coup_s if s==1, col
		list gwf_case cow year lagmil mil  dem coup_s if e(sample) & coup_s==1 & dem==1 , clean noobs
	* post-cold war *
	qui:logit dem coup_s prev_dem colbrit ind gdppc chgdppc time time2 time3 mil if cold==0, nolog
	lincom coup_s
	est store rep3
	
	 ** How does PT sample differ from GWF autocracies: (1) partial independence year in PT but not in GWF; (2) include small pop countries. 
		*small pop countries not in GWF *
		qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 if cow~=110 & cow!=347 & cow~=522 & cow!=692 & cow!=760 & cow!=910 & cow!=950, nolog
		sum gwf_duration year if e(sample)
		lincom coup_s
		* add military *
		qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 mil if  cow~=110 & cow!=347 & cow~=522 & cow!=692 & cow!=760 & cow!=910 & cow!=950, nolog
		sum gwf_duration year if e(sample)
		lincom coup_s
		* drop first years to match GWF duration>=1 *
		qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 if count~=1, nolog
		sum gwf_duration year if e(sample)
		lincom coup_s
		* add military *
		qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 mil if  count~=1, nolog
		sum gwf_duration year if e(sample)
		lincom coup_s
		* drop first years to match GWF duration>=1 & small pop countries not in GWF *
		qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 if count~=1 & cow~=110 & cow!=347 & cow~=522 & cow!=692 & cow!=760 & cow!=910 & cow!=950, nolog
		sum gwf_duration year if e(sample)
		lincom coup_s
		* add military *
		qui:logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 mil if count~=1  & cow~=110 & cow!=347 & cow~=522 & cow!=692 & cow!=760 & cow!=910 & cow!=950, nolog
		sum gwf_duration year if e(sample)
		lincom coup_s
 	

		* Show that changes to model for ANY coups also yields weaker results *
			* verification model *
			qui:logit dem coup_a prev_dem colbrit ind cold gdppc chgdppc time time2 time3, nolog
			lincom coup_a
			* add military regime *
			qui:logit dem coup_a prev_dem colbrit ind cold gdppc chgdppc time time2 time3 mil, nolog
			lincom coup_a
			* post-cold war *
			qui:logit dem coup_a prev_dem colbrit ind gdppc chgdppc time time2 time3 mil if cold==0, nolog
			lincom coup_a
 
	* linear model, PT specification *
	qui:reg dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 if gwf_duration~=. & gwf_regime~="NA", cluster(caseid)
	lincom coup_s
	est store rep4
	* linear model, add military *
	qui:reg dem mil coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 if gwf_duration~=. & gwf_regime~="NA", cluster(caseid)
	lincom coup_s
	est store rep5
	* linear model with regime-case fe
	xi: qui:reg dem i.caseid coup_s cold gdppc chgdppc time time2 time3 if gwf_duration~=. & gwf_regime~="NA", cluster(caseid)
	lincom coup_s
	est store rep6
 
	* substantive effect of adding military dummy; with post-coldwar setting *
	set seed 8979237
	qui:estsimp logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3, nolog
	setx median
	setx gdppc mean
	setx chgdppc mean
	setx cold 0  
	simqi, fd(prval(1)) changex(coup_s 0 1)
	drop b*
	qui:estsimp logit dem coup_s prev_dem colbrit ind cold gdppc chgdppc time time2 time3 mil, nolog
	setx median
	setx gdppc mean
	setx chgdppc mean
	setx cold 0  
	simqi, fd(prval(1)) changex(coup_s 0 1)
	drop b*
	qui:estsimp logit dem coup_s prev_dem colbrit ind gdppc chgdppc time time2 time3 if cold==0, nolog
	setx median
	setx gdppc mean
	setx chgdppc mean
	simqi, fd(prval(1)) changex(coup_s 0 1)
	drop b*
 
	label var coup_s "{bf:Recent coup}"
	label var gdppc "GDP pc (log)"
	label var prev_dem "Prior democracy"
	label var ind "Indep. year"
	label var colbrit "British colony"
	label var chgdppc "Growth"
	label var mil "Military regime"
	label var cold "Cold war"
	
	 * Replication plot *
	coefplot (rep1, msymbol(S)) (rep2, msymbol(T)) (rep3, msymbol(Oh)), title(Logit link, height(5))  /*
	*/ scheme(lean2) drop(_cons time* _Icase* _Iper*) xlab(-4 (1) 3) xline(0) grid(glcolor(gs15)) mfcolor(white) /*
	*/ order(coup_a prev_dem colbrit ind cold gdppc chgdppc) ysize(2) xsize(3) /*
	*/ legend(label(3 "Verification") label(6 "Add military") label(9 "post-Cold War") pos(6) ring(1.5) col(3))  /*
	*/ levels(95 90) xtitle("Coefficient estimate", height(6)) saving(t1, replace) 
	coefplot (rep4, msymbol(T)) (rep5, msymbol(Oh)) (rep6, msymbol(plus)), title("Linear link",height(5) )   /*
	*/ scheme(lean2) drop(_cons time*  _Icase*) xlab(-.04 (.01) .04) xline(0) grid(glcolor(gs15)) mfcolor(white) /*
	*/ order(coup_a prev_dem colbrit ind cold gdppc chgdppc) ysize(2) xsize(3) /*
	*/ legend(label(3 "GWF autocracies") label(6 "Add military") label(9 "Regime FE") pos(6) ring(1.5) col(3))  /*
	*/ levels(95 90) xtitle("Coefficient estimate", height(6))  byopts(xrescale) saving(t2, replace)
	graph combine t1.gph t2.gph,  xsize(7) ysize(4) graphregion(color(white))
	graph export "C:\Users\jwright\Documents\My Dropbox\Research\Coups\DFGW\Manuscript files\golden\replication.pdf", as(pdf) replace

 estout rep1 rep2 rep3 rep4 rep5 rep6 using B1.tex, cells(b(star  fmt(%9.3f)) se(par fmt(%9.2f))) stats(ll r2 N) style(tex) replace label starlevels(+ 0.10 * 0.05 ** 0.01)

*************  The End ************


log close
